Evaluating Short-Term Restoration Approaches for Degraded Grasslands in Kenya

Author

Ridhin Paul

Summary

This project evaluates the short-term effectiveness of restoration methods for degraded grasslands in Kenya. The treatments include Tilage (T), the addition of native soil inoculum (S) and Manure (M). “A grazing exclosure was installed at the campus of the Kabianga University in Kenya in April 2022, when restoration treatments were laid out and baseline soil samples taken and analysed for key soil physical, chemical and biological properties. In November 2022, fresh soil samples were taken to track short- term changes in the soil. In addition, we assessed the plant community species composition, and measured in situ soil respiration rates.” Refer “Grassland restoration experiment Kabianga Kenya.pdf” for more info.

Abrevations used

Variables and abbreviations	Explanation
M	Manure
S	Soil inoculum
T	Tillage
C	Control
YIELD	Yield in dry weight (kg/625 m^2)
DM	Dry matter (%)
Tsoil	Temperature of the topsoil (0C)
Tair	Temperature of the air (0C)
MOIS	Moisture of the topsoil (%)
pH	pH in water (dimensionless)
EC	Electrical Conductivity (µs m-1)
MBC	Microbial Biomass Carbon (mg C kg soil-1)
MBN	Microbial Biomass Nitrogen (mg N kg soil-1)
Nconc	Soil total Nitrogen (%)
Cconc	Soil total Carbon (%)
P	Olsen P (μg P g-1 soil)
BD	Bulk density (g/cm3)
NH4	Ammonium-Nitrogen (µg NH4-N g-1 DW)
NO3	Nitrate-Nitrogen (µg NO3-N g-1 DW)

Loading required packages

library(tidyverse)
library(readxl)
library(GGally)
library(VIM)
library(FactoMineR)
library(stats)
library(factoextra)
library(ggrepel)

Loading all datasets

Some unwanted columns have also been removed and -1 is replaced as 0 for easy interpretability.

#Kabianga
kab.nov.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                         sheet = 'Kab_nov_2022', 
                         na = c("NA", NA)) |>
  select(-c(...21,...22)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))
kab.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                         sheet = 'Kab_dec_2023', 
                         na = c("NA", NA)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1)) 

#Thurgem
thur.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" , 
                          sheet = 'Thur_june_2022', 
                          na = c("NA", NA)) |>
  select(-c(...13,...14)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

thur.july.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                           sheet ='Thur_july_2023', 
                           na = c("NA", NA)) |>
  select(-c(...15,...16)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

thur.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" ,
                          sheet = 'Thur_dec_2023',
                          na = c("NA", NA)) |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

#Kapsarok
kaps.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx", 
                          sheet = 'Kaps_june_2022', 
                          na = c("NA", NA)) |>
  select(-c(...13,...14))  |>
  mutate(M = if_else(M == -1, 0, 1),
         S = if_else(S == -1, 0, 1),
         T = if_else(T == -1, 0, 1))

Data exploration

Missing value check

aggr(kab.nov.22, prop = F, numbers = T)

aggr(kab.dec.23, prop = F, numbers = T)

aggr(thur.jun.22, prop = F, numbers = T)

aggr(thur.july.23, prop = F, numbers = T)

aggr(thur.dec.23, prop = F, numbers = T)

aggr(kaps.jun.22, prop = F, numbers = T)

For the dataset kab.nov.23 the plot indicate the last row has all NA entries which is removed. And for the dataset kaps.jun.22 there are 2 NA entries in the YIELD column, here the average based on the type of treatment for which the value is missing is taken and NA was replaced. Since there are only 2 NA values, code is not written for the entire column.

kab.nov.22 <- kab.nov.22|>
  na.omit()

mean_yield_1 <- kaps.jun.22 |>
  filter(Code == "S") |>
  pull(YIELD) |>
  mean(na.rm = TRUE)

mean_yield_2 <- kaps.jun.22 |>
  filter(Code == "ST") |>
  pull(YIELD) |>
  mean(na.rm = TRUE)

kaps.jun.22 <- kaps.jun.22 |>
  mutate(YIELD = if_else(Code == "S" & is.na(YIELD), mean_yield_1, YIELD),
         YIELD = if_else(Code == "ST" & is.na(YIELD), mean_yield_2, YIELD))

Checking data types

For the variables Plot No., YIELD, DM, SPE, Tsoil, Tair, MOIS, RES, pH, EC, MBC, MBN, Nconc, Cconc, P, BD, NH4, NO3, Block and for each treatment type M, S and T numeric values are expected as entry.

glimpse(kab.nov.22)

Rows: 32
Columns: 20
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil      <dbl> 18.9, 18.6, 18.7, 18.7, 18.7, 18.7, 18.8, 19.7, 20.7, 21.8,…
$ Tair       <dbl> 19.7, 19.8, 20.2, 20.3, 20.6, 20.9, 21.7, 24.6, 29.9, 31.2,…
$ MOIS       <dbl> 53.6, 37.0, 51.7, 51.6, 40.4, 52.4, 37.6, 36.6, 40.4, 38.7,…
$ YIELD      <dbl> 21.2, 13.8, 15.2, 17.8, 16.0, 16.4, 18.1, 17.1, 16.8, 15.2,…
$ NH4        <dbl> 0.7, 1.1, 2.0, 1.6, 1.2, 1.8, 0.5, 0.7, 2.1, 0.5, 1.1, 1.2,…
$ NO3        <dbl> 22.8, 31.0, 23.6, 5.3, 12.6, 24.6, 25.1, 25.3, 45.6, 14.6, …
$ P          <dbl> 1.4, 2.5, 0.6, 1.4, 1.7, 1.3, 1.0, 2.1, 1.4, 1.7, 2.0, 3.2,…
$ pH         <dbl> 5.0, 4.9, 4.7, 4.5, 4.6, 4.7, 4.7, 4.7, 4.6, 4.6, 4.4, 4.6,…
$ EC         <dbl> 41.9, 47.4, 35.1, 11.6, 19.7, 38.5, 35.2, 36.8, 59.9, 22.0,…
$ MBC        <dbl> 461.6, 512.7, 577.2, 574.4, 440.7, 514.2, 440.1, 422.5, 514…
$ MBN        <dbl> 88.0, 109.9, 108.7, 109.4, 85.9, 93.4, 81.8, 83.8, 98.5, 93…
$ Nconc      <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.4, 0.5, 0.5, 0.5, 0.4,…
$ Cconc      <dbl> 5.1, 5.4, 5.4, 5.0, 5.1, 5.3, 4.9, 4.8, 5.2, 5.3, 5.8, 4.6,…

glimpse(kab.dec.23)

Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ BD         <dbl> 0.3589024, 0.4051220, 0.3829878, 0.4227439, 0.4143293, 0.38…
$ P          <dbl> 5.9156125, 3.1895944, 1.5363629, 1.3363361, 1.8315103, 1.40…
$ YIELD      <dbl> 10.43400, 6.74075, 6.88250, 6.60400, 5.60350, 6.98075, 11.3…

glimpse(thur.jun.22)

Rows: 32
Columns: 12
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil      <dbl> 21.8, 22.5, 22.8, 23.5, 24.0, 24.4, 24.8, 25.7, 26.6, 27.7,…
$ Tair       <dbl> 24.1, 25.7, 26.7, 27.4, 28.0, 27.8, 28.7, 30.7, 32.8, 33.8,…
$ MOIS       <dbl> 44.7, 46.3, 45.0, 42.7, 28.0, 43.2, 33.3, 38.1, 33.0, 39.5,…
$ YIELD      <dbl> 4.2, 5.8, 3.8, 3.1, 3.5, 5.1, 3.9, 2.5, 3.9, 4.2, 1.6, 6.7,…
$ EC         <dbl> 136.0, 149.5, 165.3, 182.6, 176.0, 190.4, 152.0, 148.0, 128…

glimpse(thur.july.23)

Rows: 32
Columns: 14
$ Plot_No   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code      <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M         <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S         <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T         <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil     <dbl> 23.7, 23.6, 23.9, 24.8, 25.1, 25.0, 25.2, 26.7, 26.4, 28.0, …
$ Tair      <dbl> 25.3, 26.3, 28.0, 30.4, 32.3, 32.1, 32.4, 35.3, 33.5, 33.9, …
$ MOIS      <dbl> 48.1, 52.3, 47.3, 40.1, 47.9, 39.6, 47.8, 45.9, 51.6, 44.4, …
$ DW        <dbl> 162.7, 214.9, 270.3, 262.5, 298.1, 648.9, 646.1, 465.1, 253.…
$ YIELD     <dbl> 4.1, 5.4, 6.8, 6.6, 7.5, 16.2, 16.2, 11.6, 6.3, 8.4, 4.0, 19…
$ Nconc     <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, …
$ Cconc     <dbl> 1.0, 1.1, 1.3, 1.0, 1.0, 1.4, 1.2, 1.2, 0.9, 0.9, 0.8, 1.4, …

glimpse(thur.dec.23)

Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment  <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code       <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M          <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S          <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T          <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ MOIS       <dbl> 55.30603, 63.17742, 60.24427, 56.02128, 59.56154, 60.94444,…
$ P          <dbl> 2.7800613, 1.6521765, 1.6484200, 1.5918039, 1.6554478, 1.41…
$ YIELD      <dbl> 2.59225, 2.85375, 1.30200, 1.55025, 2.62850, 4.39375, 4.625…

glimpse(kaps.jun.22)

Rows: 32
Columns: 12
$ Plot_No.  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code      <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M         <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S         <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T         <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil     <dbl> 21.9, 21.9, 22.4, 23.5, 24.2, 24.7, 25.4, 26.2, 26.7, 28.7, …
$ Tair      <dbl> 24.3, 24.9, 26.1, 27.0, 28.8, 30.9, 32.5, 31.4, 32.0, 35.6, …
$ MOIS      <dbl> 12.0, 15.1, 19.9, 19.3, 25.4, 28.8, 26.2, 22.4, 18.9, 24.3, …
$ YIELD     <dbl> 1.866667, 2.900000, 1.400000, 2.300000, 10.300000, 6.500000,…
$ EC        <dbl> 106.0, 136.6, 140.8, 156.4, 110.8, 127.8, 101.8, 184.4, 140.…

All results were as expected.

Check for multicolinearity

Kabianga November 2022
```
kab.nov.22 |>
  select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
  ggpairs()
```
For the dataset a high correlation is observed between,

Var 1 Var 2 r

Tair Tsoil 0.909

Nconc Tair 0.615

YIELD MOIS 0.571

MBN NH4 0.651

EC NO3 0.980

Nconc NO3 0.603

MBN MBC 0.891

Cconc MBC 0.627

Cconc Nconc 0.731

Var 1	Var 2	r
Tair	Tsoil	0.909
Nconc	Tair	0.615
YIELD	MOIS	0.571
MBN	NH4	0.651
EC	NO3	0.980
Nconc	NO3	0.603
MBN	MBC	0.891
Cconc	MBC	0.627
Cconc	Nconc	0.731

Kabianga December 2023

kab.dec.23 |>
  select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
  ggpairs()

No significant correlations were seen for the dataset.

Thurgem June 2022
```
thur.jun.22 |>
  select(-`Plot No.`, - Treatment, -M, -S, -T) |> 
  ggpairs()
```
For the dataset a high correlation is observed between,

Var 1 Var 2 r

Tair Tsoil 0.974

Var 1	Var 2	r
Tair	Tsoil	0.974

Thurgem July 2023
```
thur.july.23 |>
  select(-Plot_No, - Treatment, -M, -S, -T) |> 
  ggpairs()
```
For the dataset a high correlation is observed between,

Var 1 Var 2 r

Tair Tsoil 0.912

YIELD DW 1

Var 1	Var 2	r
Tair	Tsoil	0.912
YIELD	DW	1

Thurgem Dec 2023

thur.dec.23 |>
  select(- `Plot No.`, - Treatment, -M, -S, -T) |> 
  ggpairs()

No significant correlations were observed between variables.

Kapsorok June 2022
```
kaps.jun.22 |>
  select(-Plot_No., - Treatment, -M, -S, -T) |> 
  ggpairs()
```
For the dataset a high correlation is observed between,

Var 1 Var 2 r

Tair Tsoil 0.962

Var 1	Var 2	r
Tair	Tsoil	0.962

Correlation between Tair and Tsoil is consistent across all datasets.

Histogram of YIELD by Treatment faceted by Block

Kabianga November 2022

ggplot(kab.nov.22, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Nov 2022)", x
       = "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

Kabianga December 2023

ggplot(kab.dec.23, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Dec 2023)", x
       = "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

Thurgem June 2022

ggplot(thur.jun.22, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem June 2022)", x
       = "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

Thurgem July 2023

ggplot(thur.july.23, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem July 2023)", x
       = "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

Thurgem Dec 2023

ggplot(thur.dec.23, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem Dec 2023)", x =
         "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

Kapsorok June 2022

ggplot(kaps.jun.22, aes(x = Code, y = YIELD)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of YIELD by Treaetment faceted by Block (Kapsorok June 2022)", x = "Treatment", y = "Yield") +
    facet_wrap(~Block) +
  theme_minimal()

For all dataset the variation of YIELD across treatment across each block seem to vary slightly however it’s statistical significance need to be checked for each dataset.

Average yield across treatment

Kabianga November 2022

kab.nov.22 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean Yield") +
  theme_minimal()

For the dataset, it can be observed that for treatments MST, MT, ST and T there is a reduction in YIELD compared to C (Control).

Kabianga December 2023

kab.dec.23 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga December 2023)", x = "Treatment", y = "Mean
       Yield") +
  theme_minimal()

For the dataset, it can be observed that for treatments MSTand M there is a reduction in YIELD compared to C (Control) and T the increase is fairly low.

Thurgem June 2022

thur.jun.22 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem June 2022)", x = "Treatment", y = "Mean Yield") +
  theme_minimal()

For the dataset, it can be observed that for treatments M, MS, MST, S, MT, ST and T there is a reduction in YIELD compared to C (Control).

Thurgem July 2023

thur.july.23 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem July 2023)", x = "Treatment", y = "Mean Yield") +
  theme_minimal()

For the dataset, it can be observed that for treatments M, MST, and ST there is a reduction in YIELD compared to C (Control) and for T a significant increase can be seen.

Thurgem Dec 2023

thur.dec.23 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem Dec 2023)", x = "Treatment", y = "Mean Yield") +
  theme_minimal()

For the dataset, it can be observed that for every treatment there is an increase in YIELD compared to treatment.

Kapsorok June 2022

kaps.jun.22 |>
  summarise(m_y_t = mean(YIELD), .by = Code) |>
  ggplot(aes(x = Code, y = m_y_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean YIELD by Treaetment (Kapsorok June 2022)", x = "Treatment", y = "Mean Yield") +
  theme_minimal()

For the dataset, it can be observed that for treatments M, MST and S there is a reduction in YIELD compared to C (Control) and for the rest of the treatments there is a significant increase in YIELD.

Visually tillage seem to have a negative effect on YIELD and manure combined with soil inoculum seems to have a slight positive effect on productivity but cannot be concluded.

Average MBC and MBN across treatment

Kabianga November 2023

It seem that manure slightly increase both MBC and MBN which is expected because MBC and MBN had a high corelation. It is also observed that tillage has a slightly negative effect.

kab.nov.22 |>
  summarise(m_mbc_t = mean(MBC), .by = Code) |>
  ggplot(aes(x = Code, y = m_mbc_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean MBC by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBC") +
  theme_minimal()

kab.nov.22 |>
  summarise(m_mbn_t = mean(MBN), .by = Code) |>
  ggplot(aes(x = Code, y = m_mbn_t)) +
  geom_bar(color = "black", alpha = 0.7, stat = "identity") +
  labs(title = "Bar chart of mean MBN by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBN") +
  theme_minimal()

Primary analysis

Common function for brute force algorithm based on AIC. This was not implemented but only tested except for the dataset Kabianga November 2022 with dependent as YIELD because it yielded a sligtly better model. Hybrid stepwise algorithm is used based on AIC keeping the interaction effects as fixed for rest of the models. And for all models p<0.10 is considered significant.

#Pass the base model as the formulae and scope as a vector of variable names under consideration

bestModel <- function(base_model, scope, dataset)
  {
    all_combs <- list()
    #seq_along creates n sequence of elements from the inputed vector
    for (i in seq_along(scope)) 
      {
          #combn creates all possible n combinations, simplify is set to false to                 obtain the result as list
          combs <- combn(scope, i, simplify = FALSE)
          
          all_combs <- c(all_combs, combs)
      }
    
    combine <- function(all_combs) 
      {
        #deparse convert type formula to string
        paste(deparse(base_model), paste(all_combs, collapse =
                                   " + "), sep = " + ")
      }
    
    #applies the function to list all_combs and returns a list which is converted to vector
    formulas <- unlist(lapply(all_combs, combine))
    
    fitted_model <- list()
    
    for (formula in formulas)
      {
        model <- lm(as.formula(formula), data = dataset)
        
        #adjust below to consider other criteria apart from AIC and also ajust the following           logic accordingly after iteration
        fitted_model[[formula]] <- list(models = model, AIC = AIC(model))
  
      }
    
    #takes in list, applies a function and returns vector
    aic_values <- sapply(fitted_model, function(x) x$AIC)
    

    #returns the fitted model and corresponding AIC
    return(fitted_model[[which.min(aic_values)]])
  
  }

Kabianga November 2022

YIELD as the dependent variable.

im <- lm(YIELD ~ T * S * M, data = kab.nov.22)

model_kab_Nov22_YIELD <- step(im, 
                              scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 +
                        P + pH + EC + MBC + MBN + Nconc + Cconc 
                    ),
                    direction = "both")

Start:  AIC=78.57
YIELD ~ T * S * M

        Df Sum of Sq    RSS    AIC
+ MOIS   1   25.6598 200.42 76.710
+ MBN    1   15.9946 210.09 78.217
<none>               226.08 78.565
+ NH4    1   11.2688 214.81 78.929
+ pH     1   10.2982 215.78 79.073
+ MBC    1    9.2920 216.79 79.222
+ Tsoil  1    7.6671 218.41 79.461
+ Tair   1    7.2481 218.83 79.523
+ EC     1    6.5224 219.56 79.628
+ NO3    1    5.6836 220.40 79.751
+ P      1    2.5974 223.49 80.195
+ Cconc  1    1.7285 224.35 80.320
+ Nconc  1    0.0001 226.08 80.565

Step:  AIC=76.71
YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
+ pH     1   13.4101 187.01 76.494
+ MBN    1   13.0343 187.39 76.558
<none>               200.42 76.710
+ P      1   10.0996 190.32 77.056
+ Tair   1    9.3024 191.12 77.189
+ NH4    1    9.1969 191.23 77.207
+ MBC    1    6.5987 193.82 77.639
+ Nconc  1    4.2286 196.19 78.028
+ Tsoil  1    3.4636 196.96 78.152
- MOIS   1   25.6598 226.08 78.565
+ Cconc  1    0.0500 200.37 78.702
+ NO3    1    0.0310 200.39 78.705
+ EC     1    0.0071 200.42 78.709

Step:  AIC=76.49
YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
<none>               187.01 76.494
+ Tair   1   11.2434 175.77 76.510
- pH     1   13.4101 200.42 76.710
+ Tsoil  1    8.3602 178.65 77.031
+ MBN    1    7.6768 179.34 77.153
+ NH4    1    5.8166 181.20 77.483
+ Nconc  1    4.1696 182.84 77.773
+ P      1    4.1162 182.90 77.782
+ MBC    1    2.3911 184.62 78.082
+ Cconc  1    0.0688 186.94 78.482
+ EC     1    0.0285 186.98 78.489
+ NO3    1    0.0031 187.01 78.494
- MOIS   1   28.7718 215.78 79.073

summary(model_kab_Nov22_YIELD)


Call:
lm(formula = YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + 
    T:S:M, data = kab.nov.22)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9018 -1.4359  0.3215  1.4152  7.5303 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  18.6138    14.5300   1.281   0.2135  
T             1.2088     2.8672   0.422   0.6774  
S             5.0615     2.2751   2.225   0.0367 *
M             1.7076     2.0662   0.826   0.4174  
MOIS          0.2926     0.1590   1.840   0.0793 .
pH           -3.4092     2.7143  -1.256   0.2223  
T:S          -3.2100     3.0478  -1.053   0.3037  
T:M          -2.1374     2.9200  -0.732   0.4719  
S:M          -6.6686     3.0928  -2.156   0.0423 *
T:S:M         5.2884     4.3188   1.225   0.2337  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.916 on 22 degrees of freedom
Multiple R-squared:  0.4955,    Adjusted R-squared:  0.2891 
F-statistic: 2.401 on 9 and 22 DF,  p-value: 0.04534

#implementation of brute force algorithm
scope = c("Tsoil", "Tair", "MOIS", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc")

base_model = YIELD ~ T * S * M

#run only if necessar of running bcz it takes a quite a bit of time the model m below is the result of this function run.
# m <- bestModel(base_model, scope, kab.nov.22)

# summary(m$model)

#model m below is the result of this function run (best_Model)
m <- lm(YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22) 
summary(m)


Call:
lm(formula = YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3236 -1.3683 -0.1658  1.0606  4.9609 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -1.38995    8.65073  -0.161  0.87396   
T             3.71723    2.81815   1.319  0.20206   
S             6.29647    2.10707   2.988  0.00727 **
M             1.82889    1.95230   0.937  0.36004   
MOIS          0.42629    0.16314   2.613  0.01665 * 
NH4           1.79159    1.49889   1.195  0.24596   
MBN           0.12320    0.06815   1.808  0.08572 . 
Nconc       -37.25304   15.70888  -2.371  0.02787 * 
T:S          -2.37206    2.73333  -0.868  0.39578   
T:M          -2.22589    2.72349  -0.817  0.42339   
S:M          -7.41175    2.84954  -2.601  0.01709 * 
T:S:M         3.27673    3.82547   0.857  0.40184   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.695 on 20 degrees of freedom
Multiple R-squared:  0.6081,    Adjusted R-squared:  0.3926 
F-statistic: 2.821 on 11 and 20 DF,  p-value: 0.0212

These results indicate that soil inoculum and the soil incolum combined with manure, moisture, MBN, and nitrogen concentration significantly influence yield, while other factors and interactions have less impact. Soil inoculum had a significant positive impact while soil incolum combined with manure had a negative impact.

MBC as the dependent variable

im <- lm(MBC ~ T * S * M, data = kab.nov.22)

model_kab_Nov22_MBC <- step(im, 
                              scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 +
                        P + pH + EC  + MBN + Nconc + Cconc 
                    ),
                    direction = "both")

Start:  AIC=261.84
MBC ~ T * S * M

        Df Sum of Sq   RSS    AIC
+ MBN    1     55169 14274 213.21
+ Nconc  1     26477 42966 248.48
+ Cconc  1     23970 45473 250.29
+ NH4    1     16758 52685 255.00
+ EC     1     14522 54921 256.33
+ Tair   1     13595 55848 256.87
+ Tsoil  1     13523 55920 256.91
+ NO3    1     13448 55995 256.95
+ pH     1      5741 63702 261.08
<none>               69443 261.84
+ P      1      3725 65718 262.08
+ MOIS   1       654 68789 263.54

Step:  AIC=213.21
MBC ~ T + S + M + MBN + T:S + T:M + S:M + T:S:M

        Df Sum of Sq   RSS    AIC
+ P      1      2644 11629 208.66
+ Cconc  1      1248 13026 212.29
<none>               14274 213.21
+ NO3    1       819 13455 213.32
+ EC     1       533 13740 214.00
+ Nconc  1       519 13754 214.03
+ pH     1       349 13925 214.42
+ Tair   1       234 14040 214.69
+ NH4    1       119 14154 214.94
+ Tsoil  1        75 14199 215.05
+ MOIS   1        49 14224 215.10
- MBN    1     55169 69443 261.84

Step:  AIC=208.66
MBC ~ T + S + M + MBN + P + T:S + T:M + S:M + T:S:M

        Df Sum of Sq   RSS    AIC
+ pH     1      1704  9925 205.59
<none>               11629 208.66
+ NH4    1       475 11154 209.32
+ Cconc  1       312 11318 209.79
+ Tair   1        68 11561 210.47
+ MOIS   1        60 11569 210.49
+ Nconc  1        22 11607 210.60
+ Tsoil  1         7 11622 210.64
+ EC     1         4 11625 210.65
+ NO3    1         0 11629 210.66
- P      1      2644 14274 213.21
- MBN    1     54089 65718 262.08

Step:  AIC=205.59
MBC ~ T + S + M + MBN + P + pH + T:S + T:M + S:M + T:S:M

        Df Sum of Sq   RSS    AIC
+ NH4    1       741  9184 205.10
<none>                9925 205.59
+ Cconc  1       382  9543 206.33
+ Tsoil  1       231  9694 206.83
+ Tair   1       219  9706 206.87
+ MOIS   1        58  9867 207.40
+ Nconc  1         5  9920 207.57
+ EC     1         3  9922 207.58
+ NO3    1         1  9924 207.58
- pH     1      1704 11629 208.66
- P      1      4000 13925 214.42
- MBN    1     44873 54798 258.26

Step:  AIC=205.1
MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + S:M + T:S:M

        Df Sum of Sq   RSS    AIC
<none>                9184 205.10
+ Cconc  1       487  8697 205.36
- NH4    1       741  9925 205.59
+ Tair   1       225  8959 206.31
+ Tsoil  1       160  9024 206.54
+ MOIS   1        74  9111 206.85
+ Nconc  1        26  9159 207.01
+ NO3    1        14  9170 207.06
+ EC     1        10  9174 207.07
- pH     1      1970 11154 209.32
- P      1      4606 13790 216.11
- MBN    1     36389 45573 254.36

summary(model_kab_Nov22_MBC)


Call:
lm(formula = MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + 
    S:M + T:S:M, data = kab.nov.22)

Residuals:
   Min     1Q Median     3Q    Max 
-29.70 -10.73  -3.20  10.10  40.76 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 368.5011   126.1290   2.922  0.00844 ** 
T           -21.6695    16.0188  -1.353  0.19123    
S            10.4362    16.8553   0.619  0.54280    
M            -4.2831    15.4943  -0.276  0.78505    
MBN           4.3856     0.4927   8.902 2.15e-08 ***
P           -22.8369     7.2111  -3.167  0.00485 ** 
pH          -46.2661    22.3393  -2.071  0.05151 .  
NH4         -14.6657    11.5453  -1.270  0.21856    
T:S         -36.6072    24.2079  -1.512  0.14612    
T:M           2.4122    21.8744   0.110  0.91329    
S:M         -17.4707    23.1437  -0.755  0.45912    
T:S:M        41.7015    33.6025   1.241  0.22896    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.43 on 20 degrees of freedom
Multiple R-squared:  0.9063,    Adjusted R-squared:  0.8548 
F-statistic: 17.59 on 11 and 20 DF,  p-value: 6.145e-08

MBN as the dependent variable

im <- lm(MBN ~ T * S * M, data = kab.nov.22)

model_kab_Nov22_MBN <- step(im, 
                              scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + Block + Tsoil + Tair + MOIS + NH4 + NO3 +
                        P + pH + EC  + MBC + Nconc + Cconc 
                    ),
                    direction = "both")

Start:  AIC=161.23
MBN ~ T * S * M

        Df Sum of Sq     RSS    AIC
+ MBC    1   2378.00  615.24 112.60
+ Nconc  1   1137.18 1856.06 147.94
+ Block  1   1051.14 1942.10 149.38
+ NH4    1   1037.41 1955.83 149.61
+ Cconc  1    846.47 2146.77 152.59
+ Tsoil  1    640.22 2353.02 155.53
+ Tair   1    573.91 2419.34 156.42
+ EC     1    538.12 2455.12 156.89
+ NO3    1    434.67 2558.57 158.21
<none>               2993.24 161.23
+ pH     1    180.42 2812.82 161.24
+ MOIS   1     18.69 2974.55 163.03
+ P      1      5.05 2988.19 163.17

Step:  AIC=112.6
MBN ~ T + S + M + MBC + T:S + T:M + S:M + T:S:M

        Df Sum of Sq     RSS    AIC
+ NH4    1     89.79  525.45 109.55
+ P      1     86.46  528.79 109.75
<none>                615.24 112.60
+ Block  1     28.78  586.47 113.07
+ Nconc  1     21.08  594.17 113.49
+ Tsoil  1     17.77  597.47 113.66
+ Tair   1      7.04  608.20 114.23
+ EC     1      1.02  614.23 114.55
+ NO3    1      0.46  614.78 114.58
+ pH     1      0.38  614.87 114.58
+ Cconc  1      0.30  614.94 114.58
+ MOIS   1      0.17  615.08 114.59
- MBC    1   2378.00 2993.24 161.23

Step:  AIC=109.55
MBN ~ T + S + M + MBC + NH4 + T:S + T:M + S:M + T:S:M

        Df Sum of Sq     RSS    AIC
+ P      1    103.39  422.07 104.54
<none>                525.45 109.55
+ NO3    1     10.27  515.19 110.92
+ Block  1      8.73  516.72 111.02
+ Tsoil  1      4.64  520.81 111.27
+ EC     1      2.46  522.99 111.40
+ Tair   1      1.68  523.78 111.45
+ Cconc  1      1.40  524.05 111.47
+ Nconc  1      1.13  524.32 111.48
+ pH     1      0.97  524.49 111.49
+ MOIS   1      0.36  525.09 111.53
- NH4    1     89.79  615.24 112.60
- MBC    1   1430.38 1955.83 149.61

Step:  AIC=104.54
MBN ~ T + S + M + MBC + NH4 + P + T:S + T:M + S:M + T:S:M

        Df Sum of Sq     RSS    AIC
+ pH     1     40.79  381.28 103.29
+ Nconc  1     28.27  393.80 104.32
+ Tair   1     26.83  395.24 104.44
<none>                422.07 104.54
+ Block  1     19.79  402.27 105.00
+ Tsoil  1     10.84  411.23 105.71
+ EC     1      7.14  414.92 106.00
+ MOIS   1      4.89  417.17 106.17
+ NO3    1      3.66  418.40 106.26
+ Cconc  1      1.65  420.41 106.42
- P      1    103.39  525.45 109.55
- NH4    1    106.72  528.79 109.75
- MBC    1   1518.73 1940.80 151.36

Step:  AIC=103.29
MBN ~ T + S + M + MBC + NH4 + P + pH + T:S + T:M + S:M + T:S:M

        Df Sum of Sq     RSS    AIC
+ Tair   1     35.09  346.19 102.20
+ Tsoil  1     24.36  356.91 103.18
+ Block  1     23.93  357.35 103.22
<none>                381.28 103.29
+ Nconc  1     17.90  363.38 103.75
- pH     1     40.79  422.07 104.54
+ MOIS   1      4.98  376.30 104.87
+ EC     1      4.71  376.57 104.89
+ NO3    1      2.85  378.43 105.05
+ Cconc  1      0.06  381.21 105.28
- NH4    1    117.87  499.14 109.91
- P      1    143.21  524.49 111.49
- MBC    1   1510.65 1891.93 152.55

Step:  AIC=102.2
MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + T:M + S:M + 
    T:S:M

        Df Sum of Sq     RSS    AIC
<none>                346.19 102.20
- Tair   1     35.09  381.28 103.29
+ MOIS   1      7.22  338.97 103.53
+ Nconc  1      3.01  343.18 103.92
+ Tsoil  1      2.08  344.11 104.01
+ Block  1      1.84  344.35 104.03
+ NO3    1      1.74  344.45 104.04
+ EC     1      1.38  344.81 104.07
+ Cconc  1      0.21  345.98 104.18
- pH     1     49.05  395.24 104.44
- NH4    1    102.43  448.62 108.49
- P      1    176.78  522.97 113.40
- MBC    1   1297.49 1643.68 150.05

summary(model_kab_Nov22_MBN)


Call:
lm(formula = MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + 
    T:M + S:M + T:S:M, data = kab.nov.22)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.6871 -2.5162  0.5216  2.5615  5.6404 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.0872    28.6689  -1.363   0.1887    
T             2.3904     3.3107   0.722   0.4791    
S            -2.6106     3.4027  -0.767   0.4524    
M             2.5024     3.0805   0.812   0.4267    
MBC           0.1746     0.0207   8.439 7.52e-08 ***
NH4           4.9858     2.1028   2.371   0.0285 *  
P             5.0085     1.6079   3.115   0.0057 ** 
pH            7.6942     4.6897   1.641   0.1173    
Tair         -0.2306     0.1661  -1.388   0.1813    
T:S           8.2018     4.9141   1.669   0.1115    
T:M          -1.4562     4.3471  -0.335   0.7413    
S:M           2.1094     4.6662   0.452   0.6563    
T:S:M        -7.2801     6.8728  -1.059   0.3028    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.269 on 19 degrees of freedom
Multiple R-squared:  0.909, Adjusted R-squared:  0.8515 
F-statistic: 15.81 on 12 and 19 DF,  p-value: 2.073e-07

MBC and MBN did not show any significant effects on treatment, but both had a very high significant effect on Phosphorus . MBC had a negative effect, whereas MBN had a positive effect on P.

Kabianga December 2023

im <- lm(YIELD ~ T * S * M, data = kab.dec.23)

model_kab_Dec23_YIELD <- step(im, 
                              scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + P + BD 
                    ),
                    direction = "both")

Start:  AIC=59.16
YIELD ~ T * S * M

       Df Sum of Sq    RSS    AIC
<none>              123.28 59.159
+ P     1   2.30009 120.98 60.556
+ BD    1   0.36336 122.92 61.064

summary(model_kab_Dec23_YIELD)


Call:
lm(formula = YIELD ~ T * S * M, data = kab.dec.23)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6739 -1.2265 -0.0384  1.1020  3.6000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.1330     1.1332   7.177 2.04e-07 ***
T             0.6874     1.6026   0.429   0.6718    
S             3.1477     1.6026   1.964   0.0612 .  
M            -0.7659     1.6026  -0.478   0.6370    
T:S          -1.2792     2.2664  -0.564   0.5777    
T:M           1.2228     2.2664   0.540   0.5945    
S:M          -1.2809     2.2664  -0.565   0.5772    
T:S:M        -2.4089     3.2052  -0.752   0.4596    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.266 on 24 degrees of freedom
Multiple R-squared:  0.3125,    Adjusted R-squared:  0.112 
F-statistic: 1.559 on 7 and 24 DF,  p-value: 0.1958

Soil inoculum showed a positive effect, which was significant and soil inoculum combined with manure had a negative effect, even though not statistically significant. The effect of treatment seems comparable to the initial measurement, even though a reduction in YIELD is observed.

Thurgem June 2022

im <- lm(YIELD ~ T * S * M, data = thur.jun.22)

model_Thur_Jun22 <- step(im,
                        scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + Tair + Tsoil + MOIS + EC 
                    ),
                    direction = "both")

Start:  AIC=26.79
YIELD ~ T * S * M

        Df Sum of Sq    RSS    AIC
+ MOIS   1   11.4629 33.370 19.341
<none>               44.832 26.790
+ Tair   1    0.2215 44.611 28.632
+ EC     1    0.1946 44.638 28.651
+ Tsoil  1    0.1631 44.669 28.674

Step:  AIC=19.34
YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
<none>               33.370 19.341
+ EC     1    0.4063 32.963 20.949
+ Tair   1    0.1978 33.172 21.151
+ Tsoil  1    0.0801 33.289 21.264
- MOIS   1   11.4629 44.832 26.790

summary(model_Thur_Jun22)


Call:
lm(formula = YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M, 
    data = thur.jun.22)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4725 -0.3853 -0.0512  0.5106  2.5927 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.25689    1.23111   1.021  0.31790   
T           -0.70152    0.85285  -0.823  0.41920   
S           -0.28573    0.85199  -0.335  0.74039   
M           -1.33047    0.85187  -1.562  0.13198   
MOIS         0.08097    0.02881   2.811  0.00992 **
T:S          0.90545    1.21536   0.745  0.46382   
T:M          3.16417    1.20452   2.627  0.01507 * 
S:M          1.60516    1.20622   1.331  0.19632   
T:S:M       -3.33147    1.70502  -1.954  0.06298 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.205 on 23 degrees of freedom
Multiple R-squared:  0.4221,    Adjusted R-squared:  0.2211 
F-statistic:   2.1 on 8 and 23 DF,  p-value: 0.07843

Tillage, when combined with manure, showed a significant positive effect. Moisture seems to have a very high significance with a p-value of 0.00992 with a positive effect.

Thurgem July 2023

Variable DW was avoided because of perfect correlation with YIELD.

im <- lm(YIELD ~ T * S * M, data = thur.july.23)

model_Thur_Jul23 <- step(im,
                        scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + Tair + Tsoil + MOIS + Nconc + Cconc +
                        Block
                    ),
                    direction = "both")

Start:  AIC=120.51
YIELD ~ T * S * M

        Df Sum of Sq    RSS    AIC
+ Cconc  1    375.54 463.02 103.50
+ Tair   1    107.83 730.74 118.11
<none>               838.56 120.51
+ MOIS   1     14.73 823.83 121.94
+ Tsoil  1     11.83 826.73 122.06
+ Block  1      1.58 836.98 122.45
+ Nconc  1      0.16 838.40 122.50

Step:  AIC=103.51
YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
<none>               463.02 103.50
+ Tair   1     25.27 437.75 103.71
+ Nconc  1     23.34 439.68 103.85
+ MOIS   1     13.93 449.09 104.53
+ Block  1      8.86 454.16 104.89
+ Tsoil  1      0.89 462.13 105.44
- Cconc  1    375.54 838.56 120.51

summary(model_Thur_Jul23)


Call:
lm(formula = YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M, 
    data = thur.july.23)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3913 -3.3870  0.1161  3.1192  7.5971 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -17.675      6.975  -2.534 0.018544 *  
T              9.336      3.417   2.732 0.011889 *  
S              3.218      3.236   0.995 0.330326    
M              6.331      3.577   1.770 0.089949 .  
Cconc         21.942      5.080   4.319 0.000254 ***
T:S           -9.041      4.551  -1.987 0.059001 .  
T:M           -9.553      4.874  -1.960 0.062238 .  
S:M           -5.955      4.826  -1.234 0.229726    
T:S:M          5.002      6.663   0.751 0.460427    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.487 on 23 degrees of freedom
Multiple R-squared:  0.5609,    Adjusted R-squared:  0.4082 
F-statistic: 3.673 on 8 and 23 DF,  p-value: 0.006751

The treatment tillage had a significant positive effect, whereas tillage combined with soil inoculum and manure had a significant negative effect. Manure also showed a significant positive effect, although its effect is not as pronounced as tillage, with an estimate of 6.331, whereas tillage had 9.336. Carbon concentration had a very significance with a p-value of 0.000254 positive effect.

Thurgem Dec 2023

im <- lm(YIELD ~ T * S * M , data = thur.dec.23)

model_Thur_dec23 <- step(im,
                        scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + MOIS + P
                    ),
                    direction = "both")

Start:  AIC=10.88
YIELD ~ T * S * M

       Df Sum of Sq    RSS    AIC
<none>              27.267 10.878
+ MOIS  1   0.64431 26.622 12.113
+ P     1   0.46632 26.800 12.326

summary(model_Thur_dec23)


Call:
lm(formula = YIELD ~ T * S * M, data = thur.dec.23)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6688 -0.6007 -0.1421  0.6171  2.6133 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.4113     0.5329   4.524 0.000139 ***
T             1.0717     0.7537   1.422 0.167918    
S             0.6397     0.7537   0.849 0.404415    
M             0.8078     0.7537   1.072 0.294468    
T:S          -0.5609     1.0659  -0.526 0.603575    
T:M          -0.7722     1.0659  -0.724 0.475789    
S:M           0.2662     1.0659   0.250 0.804918    
T:S:M        -0.5409     1.5074  -0.359 0.722870    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.066 on 24 degrees of freedom
Multiple R-squared:  0.1976,    Adjusted R-squared:  -0.03645 
F-statistic: 0.8443 on 7 and 24 DF,  p-value: 0.5624

Measurement done after a month (July 2023) showed no significant effects for any of the treatments, and the model obtained was fairly poor, with an R-squared of 19.76%.

Kapsorok June 2022

im <- lm(YIELD ~ T * S * M , data = kaps.jun.22)

model_Kaps_jun22 <- step(im,
                        scope = list(
                      lower = ~ T * S * M, 
                      upper = ~ T * S * M + MOIS + Tsoil + Tair + EC + Block
                    ),
                    direction = "both")

Start:  AIC=51.3
YIELD ~ T * S * M

        Df Sum of Sq    RSS    AIC
+ Block  1   20.3776 76.072 45.710
<none>               96.449 51.305
+ EC     1    5.3695 91.080 51.472
+ MOIS   1    0.8644 95.585 53.017
+ Tair   1    0.5081 95.941 53.136
+ Tsoil  1    0.0257 96.423 53.296

Step:  AIC=45.71
YIELD ~ T + S + M + Block + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
+ MOIS   1   13.2794 62.792 41.571
<none>               76.072 45.710
+ Tair   1    0.1182 75.953 47.660
+ Tsoil  1    0.0763 75.995 47.678
+ EC     1    0.0501 76.022 47.689
- Block  1   20.3776 96.449 51.305

Step:  AIC=41.57
YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + T:S:M

        Df Sum of Sq    RSS    AIC
<none>               62.792 41.571
+ EC     1     0.585 62.207 43.271
+ Tair   1     0.336 62.456 43.399
+ Tsoil  1     0.242 62.551 43.448
- MOIS   1    13.279 76.072 45.710
- Block  1    32.793 95.585 53.017

summary(model_Kaps_jun22)


Call:
lm(formula = YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + 
    T:S:M, data = kaps.jun.22)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9130 -0.5417 -0.0770  0.4879  4.2881 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.88600    1.34779   2.141  0.04358 * 
T            2.15809    1.19620   1.804  0.08492 . 
S           -1.10692    1.19598  -0.926  0.36473   
M           -0.95527    1.19678  -0.798  0.43329   
Block       -1.35670    0.40025  -3.390  0.00264 **
MOIS         0.12675    0.05876   2.157  0.04220 * 
T:S          1.18085    1.69415   0.697  0.49309   
T:M          0.06027    1.73371   0.035  0.97258   
S:M          3.45978    1.71023   2.023  0.05539 . 
T:S:M       -4.98098    2.49194  -1.999  0.05813 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.689 on 22 degrees of freedom
Multiple R-squared:  0.594, Adjusted R-squared:  0.4279 
F-statistic: 3.577 on 9 and 22 DF,  p-value: 0.007063

Measurements showed a significant negative effect for tillage combined with soil inoculum and manure, tillage showed a positive effect, and soi inoculum combined with manure showed a positive effect. Moisture showed a positive effect with a p-value of 0.04220.

Testing the significance of Block on dependent variable

summary(aov(YIELD ~ Block, data = kab.nov.22))

            Df Sum Sq Mean Sq F value Pr(>F)
Block        1    9.0   8.978   0.745  0.395
Residuals   30  361.7  12.056

summary(aov(YIELD ~ Block, data = kab.dec.23))

            Df Sum Sq Mean Sq F value Pr(>F)
Block        1    1.8   1.801   0.304  0.585
Residuals   30  177.5   5.917

summary(aov(YIELD ~ Block, data = thur.jun.22))

            Df Sum Sq Mean Sq F value Pr(>F)
Block        1   2.58   2.576   1.401  0.246
Residuals   30  55.17   1.839

summary(aov(YIELD ~ Block, data = thur.july.23))

            Df Sum Sq Mean Sq F value Pr(>F)
Block        1    1.6    1.58   0.045  0.833
Residuals   30 1053.0   35.10

summary(aov(YIELD ~ Block, data = thur.dec.23))

            Df Sum Sq Mean Sq F value Pr(>F)
Block        1   1.43   1.426   1.314  0.261
Residuals   30  32.56   1.085

summary(aov(YIELD ~ Block, data = kaps.jun.22))

            Df Sum Sq Mean Sq F value Pr(>F)  
Block        1  20.38  20.378   4.552 0.0412 *
Residuals   30 134.29   4.476                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Except for Kapsorok Block did not have a significant effect on YIELD. The effect of Block for Kapsorok is significant with a p-value of 0.0412.

Creating table to summarize results for report.

stargazer(model_Kaps_jun22,
          type = "text", 
          title = "Comparison of Regression Models across datssets with YIELD as dependent
          variable")

Secondary analysis

Considering the variable count the datasets on Kabianga November 2022 was only chosen for principal component analysis.

Kabianga November 2023

Performing PCA and selecting principal components

pca_result.kab <- kab.nov.22 |>
  select(-`Plot No.`, -M, -S, -T, - Treatment, -Code, -Block) |>
  PCA(scale.unit = TRUE)

fviz_eig(pca_result.kab, addlabels = TRUE, ylim = c(0, 50))

Corelation table and biplot

var_coord <- pca_result.kab$var$coord[, 1:2]

selected_vars <- var_coord[c("Tsoil", "Tair", "MOIS", "YIELD", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc"), ]

correlation_table <- data.frame(Variable = rownames(selected_vars), round(selected_vars, 2))

corelation_tibble <- as_tibble(correlation_table) |>
  rename(PC1 = Dim.1,
         PC2 = Dim.2)


print(corelation_tibble)

# A tibble: 13 × 3
   Variable   PC1   PC2
   <chr>    <dbl> <dbl>
 1 Tsoil    -0.61  0.52
 2 Tair     -0.69  0.5 
 3 MOIS      0.33  0.71
 4 YIELD     0.2   0.74
 5 NH4       0.67  0.18
 6 NO3       0.72 -0.27
 7 P        -0.49  0.22
 8 pH       -0.02 -0.06
 9 EC        0.76 -0.24
10 MBC       0.75  0.41
11 MBN       0.75  0.31
12 Nconc     0.88 -0.08
13 Cconc     0.8   0.19

K-means clustering based on principal components

Result may vary from report because set.seed(1234567) was added after report was written.

pc <- data.frame(pca_result.kab$ind$coord[, 1:2])
fviz_nbclust(pc, FUNcluster = kmeans, method = "wss")

set.seed(12334567)
k_means <- kmeans(pc, centers = 4)

fviz_pca_biplot(pca_result.kab,
                habillage = as.factor(k_means$cluster),
                addEllipses = TRUE
               )

fviz_pca_biplot(pca_result.kab,
                habillage = as.factor(kab.nov.22$Code)
               )

p_cat <- tibble(
  Code = kab.nov.22$Code,
  Cluster = as.factor(k_means$cluster)
)

ggplot(p_cat, aes(x = Code, fill=Cluster)) +
  geom_bar(stat = "count") +
  theme_minimal()

It seems as though manure treatment is associated with improved soil properties which are positively correlated with YIELD, except for EC and NO_3, which seem to be poorly correlated. Tilage seems to have a very little to no impact on YIELD.

Just a plot to visualise variables and their correlation with PC1 and PC2

correlation_data <- data.frame(
  Variable = rownames(var_coord), 
  PC1 = pca_result.kab$var$coord[, 1], 
  PC2 = pca_result.kab$var$coord[, 2]
)



# Correlation plot
ggplot(correlation_data, aes(x = PC1, y = PC2, label = Variable)) +
  geom_point(size = 3, color = "blue") +  
  geom_text(vjust = 1.5, color = "black") +  
  theme_minimal() +
  labs(title = "Correlations of Variables with Principal Components",
       x = "PC1 (Principal Component 1)",
       y = "PC2 (Principal Component 2)")

Assesing the effect of solar radiation flux and precipitation

Creating tibbble for analysis

wi_kab_nov_22 <- read_csv("data/Weather data/wi_kab_dec_23.csv")
wi_kab_dec_23 <- read_csv("data/Weather data/wi_kab_nov_22.csv")


wi_thur_jun_22 <- read_csv("data/Weather data/wi_thur_jun_22.csv")
wi_thur_jul_23 <- read_csv("data/Weather data/wi_thur_jul_23.csv")
wi_thur_dec_23 <- read_csv("data/Weather data/wi_thur_dec_23.csv")

wi_kaps_jun_22 <- read_csv("data/Weather data/wi_kaps_jun_22.csv")


a <- kab.nov.22 |>
  select(YIELD, Block, Code, MOIS) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kab_nov_22$Precipitation),
         Mean_srf = mean(wi_kab_nov_22$solar_radiation_flux),
         date = "nov_22")


b <- kab.dec.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kab_dec_23$Precipitation),
         Mean_srf = mean(wi_kab_dec_23$solar_radiation_flux),
         date = "dec_23")

c <- thur.jun.22 |>
  select(YIELD, Block, Code, MOIS) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_jun_22$Precipitation),
         Mean_srf = mean(wi_thur_jun_22$solar_radiation_flux),
         date = "jun_22")


d <- thur.july.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_jul_23$Precipitation),
         Mean_srf = mean(wi_thur_jul_23$solar_radiation_flux),
         date = "jul_23")

e <- thur.dec.23 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_thur_dec_23$Precipitation),
         Mean_srf = mean(wi_thur_dec_23$solar_radiation_flux),
         date = "dec_23")



weather_info_kab <- bind_rows(a, b)

weather_info_thur <- bind_rows(c, d ,e)

weather_info_kaps <- kaps.jun.22 |>
  select(YIELD, Block, Code) |>
  summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
  mutate(Mean_precip = mean(wi_kaps_jun_22$Precipitation),
         Mean_srf = mean(wi_kaps_jun_22$solar_radiation_flux),
         date = "jun_22")

rm(a,b,c,d,e)

Kabianga

ggplot(weather_info_kab, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) +
  geom_point() +
  geom_text_repel(aes(label = date)) +
  facet_wrap(~Code) +
  theme_minimal() +
  labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean precipitation", y = "Yield")

ggplot(weather_info_kab, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) +
  geom_line(aes(linetype = Code)) +
  geom_point() +
  geom_text_repel(aes(label = date)) +
  facet_wrap(~Code) +
  theme_minimal() +
  labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean solar radiation flux (log scale)", y = "Yield")

summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_kab))


Call:
lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_kab)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59945 -1.16581  0.05666  0.87641  2.26297 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          35.7426     2.5642  13.939 2.31e-06 ***
Mean_precip:CodeC    -5.3420     0.6415  -8.328 7.05e-05 ***
Mean_precip:CodeM    -5.2672     0.6415  -8.211 7.71e-05 ***
Mean_precip:CodeMS   -5.1600     0.6415  -8.044 8.80e-05 ***
Mean_precip:CodeMST  -5.6693     0.6415  -8.838 4.80e-05 ***
Mean_precip:CodeMT   -5.4248     0.6415  -8.457 6.38e-05 ***
Mean_precip:CodeS    -4.5853     0.6415  -7.148 0.000186 ***
Mean_precip:CodeST   -5.1048     0.6415  -7.958 9.43e-05 ***
Mean_precip:CodeT    -5.4642     0.6415  -8.518 6.09e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.769 on 7 degrees of freedom
Multiple R-squared:  0.9279,    Adjusted R-squared:  0.8454 
F-statistic: 11.25 on 8 and 7 DF,  p-value: 0.002275

summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_kab))


Call:
lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_kab)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0193 -0.7584 -0.0061  0.7237  1.8153 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -6.176e+01  7.236e+00  -8.535 6.02e-05 ***
Mean_srf:CodeC    3.680e-06  3.604e-07  10.211 1.86e-05 ***
Mean_srf:CodeM    3.711e-06  3.604e-07  10.298 1.76e-05 ***
Mean_srf:CodeMS   3.715e-06  3.604e-07  10.311 1.75e-05 ***
Mean_srf:CodeMST  3.597e-06  3.604e-07   9.982 2.17e-05 ***
Mean_srf:CodeMT   3.641e-06  3.604e-07  10.103 2.00e-05 ***
Mean_srf:CodeS    3.849e-06  3.604e-07  10.681 1.38e-05 ***
Mean_srf:CodeST   3.711e-06  3.604e-07  10.299 1.76e-05 ***
Mean_srf:CodeT    3.636e-06  3.604e-07  10.090 2.02e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.538 on 7 degrees of freedom
Multiple R-squared:  0.9454,    Adjusted R-squared:  0.8831 
F-statistic: 15.16 on 8 and 7 DF,  p-value: 0.0008942

Thurgum

ggplot(weather_info_thur, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) +
  geom_point() +
  geom_text_repel(aes(label = date)) +
  facet_wrap(~Code) +
  theme_minimal() +
  labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean
       precipitation", y = "Yield")

ggplot(weather_info_thur, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) +
  geom_line(aes(linetype = Code)) +
  geom_point() +
  geom_text_repel(aes(label = date)) +
  facet_wrap(~Code) +
  theme_minimal() +
  labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean
       solar radiation flux (log scale)", y = "Yield")

summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_thur))


Call:
lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_thur)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2480 -2.9100 -0.0426  2.0472  6.8268 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)            9.370      1.516   6.180 1.76e-05 ***
Mean_precip:CodeC     -2.885      1.716  -1.681    0.114    
Mean_precip:CodeM     -2.837      1.716  -1.653    0.119    
Mean_precip:CodeMS    -2.164      1.716  -1.261    0.227    
Mean_precip:CodeMST   -3.177      1.716  -1.851    0.084 .  
Mean_precip:CodeMT    -2.216      1.716  -1.291    0.216    
Mean_precip:CodeS     -2.639      1.716  -1.537    0.145    
Mean_precip:CodeST    -2.689      1.716  -1.567    0.138    
Mean_precip:CodeT     -2.223      1.716  -1.295    0.215    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.887 on 15 degrees of freedom
Multiple R-squared:  0.3227,    Adjusted R-squared:  -0.03853 
F-statistic: 0.8933 on 8 and 15 DF,  p-value: 0.5452

summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_thur))


Call:
lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_thur)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5003 -2.4308  0.3415  1.9385  4.8084 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       7.197e+01  1.566e+01   4.595 0.000350 ***
Mean_srf:CodeC   -3.230e-06  7.686e-07  -4.203 0.000769 ***
Mean_srf:CodeM   -3.251e-06  7.686e-07  -4.229 0.000728 ***
Mean_srf:CodeMS  -3.182e-06  7.686e-07  -4.140 0.000873 ***
Mean_srf:CodeMST -3.310e-06  7.686e-07  -4.307 0.000623 ***
Mean_srf:CodeMT  -3.169e-06  7.686e-07  -4.123 0.000903 ***
Mean_srf:CodeS   -3.215e-06  7.686e-07  -4.183 0.000799 ***
Mean_srf:CodeST  -3.242e-06  7.686e-07  -4.218 0.000745 ***
Mean_srf:CodeT   -3.165e-06  7.686e-07  -4.118 0.000913 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.093 on 15 degrees of freedom
Multiple R-squared:  0.571, Adjusted R-squared:  0.3422 
F-statistic: 2.495 on 8 and 15 DF,  p-value: 0.06039

In Kabianga, productivity was lower during the months with high rainfall. However, when considering solar radiation flux, yield increased with higher levels of solar radiation. In contrast, at Thurgem, YIELD decreased as both rainfall and solar radiation flux increased.

Summary

Abrevations used

Loading required packages

Loading all datasets

Data exploration

Missing value check

Checking data types

Check for multicolinearity

Histogram of YIELD by Treatment faceted by Block

Average yield across treatment

Average MBC and MBN across treatment

Primary analysis

Testing the significance of Block on dependent variable

Secondary analysis

Kabianga November 2023

Performing PCA and selecting principal components

Corelation table and biplot

K-means clustering based on principal components

Just a plot to visualise variables and their correlation with PC1 and PC2

Assesing the effect of solar radiation flux and precipitation

Creating tibbble for analysis

Miscellaneous